#include "mdp.h"
#include "fermiqcd_coefficients.h"

void set_cold_U1(mdp_field<mdp_complex>& U) {
  mdp_site x(U.lattice());
  forallsitesandcopies(x) 
    for(int mu=0; mu<U.lattice().ndim; mu++)
      U(x,mu)=1;
}

void set_hot_U1(mdp_field<mdp_complex>& U) {
  mdp_site x(U.lattice());
  forallsitesandcopies(x) 
    for(int mu=0; mu<U.lattice().ndim; mu++)
      U(x,mu)=exp(2.0*I*Pi*U.lattice().random(x).plain());
}

mdp_complex staple_H(mdp_field<mdp_complex>& U, 
		     mdp_site x, int mu) {
  mdp_complex tmp=0;
  for(int nu=0; nu<U.lattice().ndim; nu++) if(nu!=mu) {
    tmp+=conj(U(x,nu)*U(x+nu,mu))*U(x+mu,nu);
    tmp+=U(x-nu,nu)*conj(U(x-nu,mu))*conj(U((x+mu)-nu,nu));
  }
  return tmp;
}

mdp_complex plaquette_U1(mdp_field<mdp_complex>& U, 
			 mdp_site x, int mu, int nu) {
  return U(x,mu)*U(x+mu,nu)*conj(U(x,nu)*U(x+nu,mu));
}

mdp_real average_plaquette_U1(mdp_field<mdp_complex>& U) {
  mdp_real tmp=0;
  site x(U.lattice());
  // U.update();
  int mu,nu,ndim=U.lattice().ndim;
  forallsites(x)
    for(mu=0; mu<ndim-1; mu++)
      for(nu=mu+1; nu<ndim; nu++)
        tmp+=real(plaquette_U1(U,x,mu,nu));
  mdp.add(tmp);
  return 2.0*tmp/(ndim*(ndim-1)*U.lattice().nvol_gl);
}


void metropolis_multihit_U1(mdp_field<mdp_complex> &U, 
			 float beta, int n_iter, int n_hits) {
  int iter, parity, hit, mu, ndim=U.lattice().ndim;
  mdp_complex newU; 
  float delta_action;
  mdp_complex stapleh;
  mdp_site x(U.lattice());
  for(iter=0; iter<n_iter; iter++) {
    for(parity=0; parity<=1; parity++) {
      forallsitesofparity(x,parity) {
	for(mu=0; mu<ndim; mu++) {
	  stapleh=staple_H(U,x,mu);
	  for(hit=0; hit<n_hits; hit++) {
	    newU=exp(2.0*Pi*I*U.lattice().random(x).plain());
	    // CHECK FOLLOWING LINE FOR SIGN AND COEFFICINTS
	    delta_action=beta*real((U(x,mu)-newU)*stapleh); 
	    if(delta_action<U.lattice().random(x).plain()) U(x,mu)=newU;
	  }
	}
      }
      U.update(parity);
    }
  }
}


int main(int argc, char** argv) {
  mdp.open_wormholes(argc,argv);
  int L[]={4,4,4,4}; // the lattice size
  mdp_lattice spacetime(4,L); // in 4D
  mdp_field<mdp_complex> U(spacetime,4); 
  set_hot_U1(U);
  for(int i=0; i<10; i++) {
    metropolis_multihit_U1(U,4.0,10,1); // beta=4.0, 10 iterations, 1 hit each
    mdp << "P=" << average_plaquette_U1(U) << endl;
    // do whatever you want with it
  }
  mdp.close_wormholes();
  return 0;
}
